home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / dassl / ddajac.f < prev    next >
Text File  |  1996-07-19  |  6KB  |  178 lines
  1.       SUBROUTINE DDAJAC (NEQ, X, Y, YPRIME, DELTA, CJ, H,
  2.      +   IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR,
  3.      +   IPAR, NTEMP)
  4. C***BEGIN PROLOGUE  DDAJAC
  5. C***SUBSIDIARY
  6. C***PURPOSE  Compute the iteration matrix for DDASSL and form the
  7. C            LU-decomposition.
  8. C***LIBRARY   SLATEC (DASSL)
  9. C***TYPE      DOUBLE PRECISION (SDAJAC-S, DDAJAC-D)
  10. C***AUTHOR  PETZOLD, LINDA R., (LLNL)
  11. C***DESCRIPTION
  12. C-----------------------------------------------------------------------
  13. C     THIS ROUTINE COMPUTES THE ITERATION MATRIX
  14. C     PD=DG/DY+CJ*DG/DYPRIME (WHERE G(X,Y,YPRIME)=0).
  15. C     HERE PD IS COMPUTED BY THE USER-SUPPLIED
  16. C     ROUTINE JAC IF IWM(MTYPE) IS 1 OR 4, AND
  17. C     IT IS COMPUTED BY NUMERICAL FINITE DIFFERENCING
  18. C     IF IWM(MTYPE)IS 2 OR 5
  19. C     THE PARAMETERS HAVE THE FOLLOWING MEANINGS.
  20. C     Y        = ARRAY CONTAINING PREDICTED VALUES
  21. C     YPRIME   = ARRAY CONTAINING PREDICTED DERIVATIVES
  22. C     DELTA    = RESIDUAL EVALUATED AT (X,Y,YPRIME)
  23. C                (USED ONLY IF IWM(MTYPE)=2 OR 5)
  24. C     CJ       = SCALAR PARAMETER DEFINING ITERATION MATRIX
  25. C     H        = CURRENT STEPSIZE IN INTEGRATION
  26. C     IER      = VARIABLE WHICH IS .NE. 0
  27. C                IF ITERATION MATRIX IS SINGULAR,
  28. C                AND 0 OTHERWISE.
  29. C     WT       = VECTOR OF WEIGHTS FOR COMPUTING NORMS
  30. C     E        = WORK SPACE (TEMPORARY) OF LENGTH NEQ
  31. C     WM       = REAL WORK SPACE FOR MATRICES. ON
  32. C                OUTPUT IT CONTAINS THE LU DECOMPOSITION
  33. C                OF THE ITERATION MATRIX.
  34. C     IWM      = INTEGER WORK SPACE CONTAINING
  35. C                MATRIX INFORMATION
  36. C     RES      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
  37. C                TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME)
  38. C     IRES     = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES
  39. C                IN RES, AND LESS THAN ZERO OTHERWISE.  (IF IRES
  40. C                IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED)
  41. C                IN THIS CASE (IF IRES .LT. 0), THEN IER = 0.
  42. C     UROUND   = THE UNIT ROUNDOFF ERROR OF THE MACHINE BEING USED.
  43. C     JAC      = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE
  44. C                TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE
  45. C                IS ONLY USED IF IWM(MTYPE) IS 1 OR 4)
  46. C-----------------------------------------------------------------------
  47. C***ROUTINES CALLED  DGBFA, DGEFA
  48. C***REVISION HISTORY  (YYMMDD)
  49. C   830315  DATE WRITTEN
  50. C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
  51. C   901010  Modified three MAX calls to be all on one line.  (FNF)
  52. C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
  53. C   901026  Added explicit declarations for all variables and minor
  54. C           cosmetic changes to prologue.  (FNF)
  55. C   901101  Corrected PURPOSE.  (FNF)
  56. C***END PROLOGUE  DDAJAC
  57. C
  58.       INTEGER  NEQ, IER, IWM(*), IRES, IPAR(*), NTEMP
  59.       DOUBLE PRECISION
  60.      *   X, Y(*), YPRIME(*), DELTA(*), CJ, H, WT(*), E(*), WM(*),
  61.      *   UROUND, RPAR(*)
  62.       EXTERNAL  RES, JAC
  63. C
  64.       EXTERNAL  DGBFA, DGEFA
  65. C
  66.       INTEGER  I, I1, I2, II, IPSAVE, ISAVE, J, K, L, LENPD, LIPVT,
  67.      *   LML, LMTYPE, LMU, MBA, MBAND, MEB1, MEBAND, MSAVE, MTYPE, N,
  68.      *   NPD, NPDM1, NROW
  69.       DOUBLE PRECISION  DEL, DELINV, SQUR, YPSAVE, YSAVE
  70. C
  71.       PARAMETER (NPD=1)
  72.       PARAMETER (LML=1)
  73.       PARAMETER (LMU=2)
  74.       PARAMETER (LMTYPE=4)
  75.       PARAMETER (LIPVT=21)
  76. C
  77. C***FIRST EXECUTABLE STATEMENT  DDAJAC
  78.       IER = 0
  79.       NPDM1=NPD-1
  80.       MTYPE=IWM(LMTYPE)
  81.       GO TO (100,200,300,400,500),MTYPE
  82. C
  83. C
  84. C     DENSE USER-SUPPLIED MATRIX
  85. 100   LENPD=NEQ*NEQ
  86.       DO 110 I=1,LENPD
  87. 110      WM(NPDM1+I)=0.0D0
  88.       CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
  89.       GO TO 230
  90. C
  91. C
  92. C     DENSE FINITE-DIFFERENCE-GENERATED MATRIX
  93. 200   IRES=0
  94.       NROW=NPDM1
  95.       SQUR = SQRT(UROUND)
  96.       DO 210 I=1,NEQ
  97.          DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)),ABS(WT(I)))
  98.          DEL=SIGN(DEL,H*YPRIME(I))
  99.          DEL=(Y(I)+DEL)-Y(I)
  100.          YSAVE=Y(I)
  101.          YPSAVE=YPRIME(I)
  102.          Y(I)=Y(I)+DEL
  103.          YPRIME(I)=YPRIME(I)+CJ*DEL
  104.          CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
  105.          IF (IRES .LT. 0) RETURN
  106.          DELINV=1.0D0/DEL
  107.          DO 220 L=1,NEQ
  108. 220      WM(NROW+L)=(E(L)-DELTA(L))*DELINV
  109.       NROW=NROW+NEQ
  110.       Y(I)=YSAVE
  111.       YPRIME(I)=YPSAVE
  112. 210   CONTINUE
  113. C
  114. C
  115. C     DO DENSE-MATRIX LU DECOMPOSITION ON PD
  116. 230      CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER)
  117.       RETURN
  118. C
  119. C
  120. C     DUMMY SECTION FOR IWM(MTYPE)=3
  121. 300   RETURN
  122. C
  123. C
  124. C     BANDED USER-SUPPLIED MATRIX
  125. 400   LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ
  126.       DO 410 I=1,LENPD
  127. 410      WM(NPDM1+I)=0.0D0
  128.       CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
  129.       MEBAND=2*IWM(LML)+IWM(LMU)+1
  130.       GO TO 550
  131. C
  132. C
  133. C     BANDED FINITE-DIFFERENCE-GENERATED MATRIX
  134. 500   MBAND=IWM(LML)+IWM(LMU)+1
  135.       MBA=MIN(MBAND,NEQ)
  136.       MEBAND=MBAND+IWM(LML)
  137.       MEB1=MEBAND-1
  138.       MSAVE=(NEQ/MBAND)+1
  139.       ISAVE=NTEMP-1
  140.       IPSAVE=ISAVE+MSAVE
  141.       IRES=0
  142.       SQUR=SQRT(UROUND)
  143.       DO 540 J=1,MBA
  144.          DO 510 N=J,NEQ,MBAND
  145.           K= (N-J)/MBAND + 1
  146.           WM(ISAVE+K)=Y(N)
  147.           WM(IPSAVE+K)=YPRIME(N)
  148.           DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
  149.           DEL=SIGN(DEL,H*YPRIME(N))
  150.           DEL=(Y(N)+DEL)-Y(N)
  151.           Y(N)=Y(N)+DEL
  152. 510       YPRIME(N)=YPRIME(N)+CJ*DEL
  153.       CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
  154.       IF (IRES .LT. 0) RETURN
  155.       DO 530 N=J,NEQ,MBAND
  156.           K= (N-J)/MBAND + 1
  157.           Y(N)=WM(ISAVE+K)
  158.           YPRIME(N)=WM(IPSAVE+K)
  159.           DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N)))
  160.           DEL=SIGN(DEL,H*YPRIME(N))
  161.           DEL=(Y(N)+DEL)-Y(N)
  162.           DELINV=1.0D0/DEL
  163.           I1=MAX(1,(N-IWM(LMU)))
  164.           I2=MIN(NEQ,(N+IWM(LML)))
  165.           II=N*MEB1-IWM(LML)+NPDM1
  166.           DO 520 I=I1,I2
  167. 520         WM(II+I)=(E(I)-DELTA(I))*DELINV
  168. 530      CONTINUE
  169. 540   CONTINUE
  170. C
  171. C
  172. C     DO LU DECOMPOSITION OF BANDED PD
  173. 550   CALL DGBFA(WM(NPD),MEBAND,NEQ,
  174.      *    IWM(LML),IWM(LMU),IWM(LIPVT),IER)
  175.       RETURN
  176. C------END OF SUBROUTINE DDAJAC------
  177.       END
  178.